module compute

  use MCMC_utils, only : random_normal

  implicit none

contains

  subroutine compute_y (yi,y_current,stepsize_y,m1,z,nmeasure,&
                       dum2,Dinv,n0,dumx,dumy,Rinv_current,n0T,m0T,reject_y,&
                       reject_y_vector,timeindex_nonzero,in_subsample,burn_in_done)

    ! input only
    integer, intent(in) :: yi,nmeasure,timeindex_nonzero
    logical, intent(in) :: burn_in_done,in_subsample
    real, intent(in)    :: stepsize_y
    real, dimension(:), intent(in)    :: z
    real, dimension(:,:), intent(in)  :: Dinv,Rinv_current
    ! inout
    integer, intent(inout)            :: reject_y
    real, intent(inout)               :: n0T,m0T
    real, dimension(:), intent(inout) :: y_current,n0,reject_y_vector

    ! really locals (temporary values)
    real, dimension(:), intent(inout) :: m1,dum2,dumx,dumy

    real :: dymod, m1T,dum,dumb,dumc,dumd,n1T,pT,randomu,n1,y_new
    integer :: ii, nmeasuremax

    real, external :: sdot

    nmeasuremax = size(y_current)

    ! Generate new value of y
    dymod = random_normal()*stepsize_y
    y_new = y_current(yi) + dymod

    ! subsample at observation times
    ! so only update if yi exists in the index array
    ! must compute m1 fully and m1T in setup()
    if (in_subsample) then
      ! calculate (z-y)^T*Dinv*(z-y)
      !IB scalar update
      m1(timeindex_nonzero) = m1(timeindex_nonzero) + y_current(yi) - y_new
      !$OMP PARALLEL DO default(none), &
      !$OMP           shared(dum2, m1, Dinv, nmeasure)
      !IB vector mult (1488)
      do ii=1,nmeasure
        dum2(ii) = m1(ii)*Dinv(ii,ii)*m1(ii)
      end do
      !$OMP END PARALLEL DO
      !IB vector reduction (1488)
      m1T = sum(dum2)
    endif

    ! compute (Hx-y)^T*Rinv*(Hx-y) based on current value of y
    !
    !IB vector copy (1488)
    !n1 = n0
    !IB scalar update at yi
    n1 = n0(yi) - dymod

    !IB vector scaling (1488)
    !dumx = n0(yi)*Rinv_current(yi,:)
    call scopy (nmeasuremax, Rinv_current(1,yi), 1, dumx, 1)
    call scopy (nmeasuremax, dumx, 1, dumy, 1)
    call sscal (nmeasuremax, n0(yi), dumx, 1)
    !dumy = n1(yi)*Rinv_current(yi,:)
    call sscal (nmeasuremax, n1, dumy, 1)

    !IB vector products (1488)
    dumc = sdot(nmeasuremax,n0,1,Rinv_current(1,yi),1)
    dum  = n0(yi) * dumc
    dumd = dumc + (n1 - n0(yi)) * Rinv_current(yi,yi)
    dumb = dumd * n1
    !dum  = sdot(nmeasuremax,n0,1,dumx,1)
    !dumb = sdot(nmeasuremax,n1,1,dumy,1)
    !dumc = sdot(nmeasuremax,n0,1,Rinv_current(yi,1),nmeasuremax)
    !dumd = sdot(nmeasuremax,n1,1,Rinv_current(yi,1),nmeasuremax)
    !dum  = dot_product(n0,dumx)
    !dumb = dot_product(n1,dumy)
    !dumc = dot_product(n0,Rinv_current(yi,:))
    !dumd = dot_product(n1,Rinv_current(yi,:))

    n1T = n0T - dum + dumb - n0(yi)*dumc + n1*dumd - &
         n1*n1*Rinv_current(yi,yi) + n0(yi)*n0(yi)*Rinv_current(yi,yi)


    ! Compute P1/P0
    pT= -0.5*(m1T - m0T) -0.5*(n1T - n0T)


    call random_number(randomu)              ! Generates uniformly distributed random number

    if(alog(randomu) .le. pT) then
       !;ACCEPT
       y_current(yi) = y_new
       n0(yi)=n1
       n0T=n1T
       m0T = m1T
    else
       !;REJECT
       if(burn_in_done) then
          reject_y=reject_y+ 1
          reject_y_vector(yi)=reject_y_vector(yi)+ 1
       endif
    endif

  end subroutine compute_y
                                                                         

end module
